#### make ewic rollout map in R
#last modified: 20 November 2024
#last modified by: Charlotte Ambrozek
#this R script file produces a nicer map of the rollout timing of WIC EBT (breaks by fiscalyear)
library(haven)
library(tidyverse)
library(xtable)
library(fixest)
library(sf)
library(forcats)
library(RColorBrewer)
library(scico)

## Create regression sample summary statistics
inpath <- "./data/cleaned"
geo_inpath <- "./data/raw/tl_2014_us_county.shp"
map_outpath <- "./analysis/output/graphs"

ewic <- read_dta(file.path(inpath, "ewic_rollout.dta")) %>%
  rename(fips_code = ct_fips) %>%
  mutate(ev_year = if_else(state == "Nevada" | state == "Mississippi" | state == "Vermont", NA, ev_year)) %>% # make v year missing in VT, NV, and MS to reflect that we do not use those states in the analysis
  select(c(fips_code, ev_year)) %>%
  distinct() %>%
  group_by(fips_code)

anyDuplicated(ewic$fips_code)

county_shp <- st_read(geo_inpath) %>%
  filter(!(STATEFP %in% c("60", "66", "69", "72", "78"))) %>% # removes non-continental US (not included in analysis)
  filter(!(STATEFP %in% c("02", "15"))) %>% # removes alaska and hawaii for map purposes
  select(fips_code = GEOID, geometry) %>%
  mutate(fips_code = as.numeric(fips_code))

##  Full data availability map ----------------------------------------------------------------------
map_data <- left_join(county_shp, ewic, by = "fips_code")

ggplot(map_data) +
  geom_sf(aes(fill = ev_year), col = "black", lwd = .08) +
  theme_light() +
  labs(fill = "WIC EBT Year of First Implementation") +
  scico::scale_fill_scico(palette = "vik", breaks = seq(2002, 2022, 10)) +
  #scale_fill_distiller(palette = "RdPu") +
  guides(fill = guide_colourbar(override.aes = list(size = 1))) +
  theme(legend.position = "bottom", axis.text = element_blank(),
        axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.text = element_text(size = 8),
        plot.background = element_rect(color = "white"),
        panel.border = element_rect(color = "white"),
        legend.title = element_text(size = 9), legend.direction = "horizontal") 

ggsave(filename = file.path(map_outpath, "wicebtrollout.png"), height = 5, width = 7, dpi = 600, units = "in")